*     discharge at CZ, RTIME for areasink with corners CZ1 CZ2 CZ3 CZ4
*     in counterclockwise order with timely and spacially constant
*     extraction for RTIME>0 and 0 for RTIME<0) with transient constant
*     RP (=porosity/(hydaulic conductivity*average head)
*     it is assumed that CZ is not on one of the sides of the areasink
      COMPLEX*16 FUNCTION CFSAQ1(CZ, RTIME, CZ1, CZ2, CZ3, CZ4, RP)
      IMPLICIT NONE
      COMPLEX*16       CZ, CZ1, CZ2, CZ3, CZ4
      DOUBLE PRECISION RTIME, RP

      COMPLEX*16       CFSAQS


      IF (RTIME.LE.0.) THEN
         CFSAQ1=(0.,0.)
         RETURN
      END IF
      CFSAQ1=( CFSAQS(CZ, RTIME, CZ1, CZ2, RP)
     &        +CFSAQS(CZ, RTIME, CZ2, CZ3, RP)
     &        +CFSAQS(CZ, RTIME, CZ3, CZ4, RP)
     &        +CFSAQS(CZ, RTIME, CZ4, CZ1, RP) )/RP
      RETURN
      END



*     Qx and Qy at CZ, RTIME for line-doublet from CZ1 to CZ2 with
*     timely linear strength (1.*RTIME for RTIME>0 and 0 for
*     for RTIME<0) with transient constant RP
*     (=porosity/(hydaulic conductivity*average head))
      COMPLEX*16 FUNCTION CFSAQS(CZ, RTIME, CZ1, CZ2, RP)
      IMPLICIT NONE
      COMPLEX*16       CZ, CZ1, CZ2
      DOUBLE PRECISION RTIME, RP

      COMPLEX*16       C12, CZL
      DOUBLE PRECISION R12, RX1, RX2, RX, RY, RPO4T, RSQPO4, RLIM,
     &                 RX2XYY, RX1XYY, RQXI, RQETA

      DOUBLE PRECISION RPI, RSQRPI, RFE1, RFERFC, RFIEER

      PARAMETER (RPI=3.1415926535)
      PARAMETER (RSQRPI=1.77245385090551)


      IF (RTIME.LE.0.) THEN
         CFSAQS=0.
         RETURN
      END IF
      C12=CZ2-CZ1
      R12=ABS(C12)
*     the element is translated and rotated such that the center is at
*     the origin and the element along the x-axis with endpoint 1 at x1=
*     -.5*(length of element), endpoint 2 at x2=.5*(length of element)
      CZL=(CZ-.5*(CZ2+CZ1))*R12 / (C12)
      RX1=-.5*R12
      RX2=-RX1
      RX=REAL(CZL)
      RY=AIMAG(CZL)
      RPO4T =.25*RP/RTIME
      RSQPO4=SQRT(RPO4T)
      RLIM=RSQPO4*RY
      RX2XYY=(RX-RX2)*(RX-RX2)+RY*RY
      RX1XYY=(RX-RX1)*(RX-RX1)+RY*RY
      RQXI=.5*RY*RTIME/RPI*EXP(-RPO4T*RX2XYY)/RX2XYY
     &   -.5*RY*RTIME/RPI*EXP(-RPO4T*RX1XYY)/RX1XYY
     &   -.125*RP*RY/RPI*(RFE1(RPO4T*RX2XYY)-RFE1(RPO4T*RX1XYY))
      RQETA=-.5*RTIME/RPI*((RX-RX2)/RX2XYY*EXP(-RPO4T*RX2XYY)
     &                    -(RX-RX1)/RX1XYY*EXP(-RPO4T*RX1XYY))
     &   -.125*RP/RPI*((RX-RX2)*RFE1(RPO4T*RX2XYY)
     &                -(RX-RX1)*RFE1(RPO4T*RX1XYY))
     &   +.5*SQRT(RP*RTIME)/RSQRPI*EXP(-RPO4T*RY*RY)*
     &     (RFERFC(RSQPO4*(RX-RX2))-RFERFC(RSQPO4*(RX-RX1)))
     &   -.5*RP*RY/RSQRPI*RFIEER(RLIM, (RX-RX2)/RY, (RX-RX1)/RY)
      CFSAQS=CMPLX(RQXI, RQETA)*C12/R12
      RETURN
      END
